Numerical Study of the Localization-Delocalization Transition for Vibrations in 

Amorphous Silicon 
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Numerical studies of amorphous silicon in harmonic ap- 
proximation show that the highest 3.5% of vibrational nor- 
mal modes are localized. As vibrational frequency increases 
through the boundary separating localized from delocal- 
ized modes, near u c =70meV, (the "mobility edge") there 
is a localization-delocalization (LD) transition, similar to a 
second-order thermodynamic phase transition. By a numeri- 
cal study on a system with 4096 atoms, we are able to see ex- 
ponential decay lengths of exact vibrational eigenstates, and 
test whether or not these diverge at u c . Results are consis- 
tent with a localization length £ which diverges above lu c as 
(lj — uj c )~ p where the exponent is p ~ 1.3 ± 0.5. Below the 
mobility edge we find no evidence for a diverging correlation 
length. Such an asymmetry would contradict scaling ideas, 
and we suppose it is a finite-size artifact. If the scaling regime 
is narrower than our (~ 1 meV) resolution, then it cannot be 
seen directly on our finite system. 
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Harmonic normal modes of vibration can be classified 
as extended or localized. In dimension d = 3, the vibra- 
tional spectrum has sharp boundaries ( "mobility edges" ) 
separating these two kinds of modes. The sudden change 
at the mobility edge is a "localized to delocalized" (LD) 
transition (Anderson, 1958). Anderson's paper consid- 
ered localization of electronic wavefunctions, and most 
subsequent work focussed on this case (Lee and Ramakr- 
ishnan, 1985; Kramer and MacKinnon, 1993; Belitz and 
Kirkpatrick, 1994; Ohtsuki, Slevin, and Kawarabayashi, 
1999). Harmonic vibrations (treated classically - quanti- 
zation does not matter for localization) have similar wave 
properties to electrons in independent-electron approxi- 
mation. However, in one respect, vibrational states dif- 
fer from electron states. Long wavelength, low frequency 
(sound) modes always exist, extending throughout a ho- 
mogeneous sample, and are thus delocalized. There is no 
analog of this "gapless" mode for electrons. Neverthe- 
less, it is generally believed that at higher frequencies, 
vibrations and electrons have similar localization proper- 
ties (Kirkpatrick, 1985; Sheng, Zhou, and Zhang, 1994; 
Sheng, 1995; Akita and Ohtsuki, 1998). 



In glasses, only a few percent of normal modes can 
propagate ballistically like sound. As frequency in- 
creases, ballistic character is lost at typical frequencies of 
10 meV (Rat et ai. 1999), but the LD transition (which 
is hard to probe by direct experiment) is postponed to 
a much higher frequency uj c (Allen and Feldman 1989; 
Feldman et ai. 1991; Sheng, Zhou, and Zhang 1994; 
Finkemeier and von Niessen 1998; Pilla et ai. 2000.) Re- 
cent preprints by Galanti and Olami (1999) and Gat and 
Olami (2000) have examined the nature of vibrational 
states near the mobility edge. Using scalar vibrations, 
they find numerical evidence that vibrations may belong 
to a different "universality class" from electrons. 

The concerns of the present paper are not so much the 
formal questions of universality, but instead specific ques- 
tions of how vibrational eigenstates look in real materials. 
Experiment does not give access to vibrational eigenvec- 
tors or to electron wavefunctions. Much of the infor- 
mation about these states comes from computer studies. 
Most of these studies use simplified "toy" models, such 
as "scalar vibrations" where the vector displacement of 
an atom is replaced by a scalar variable, or hypothetical 
structures with fractal geometries (Alexander and Or- 
bach, 1982; Nakayama et ai., 1994; Kantclhardt et ai., 
1998). Here we look at properties of vibrational eigen- 
states near the mobility edge in the most realistic model 
presently available for amorphous Si. We use vector vi- 
brations and find that the localization length diverges in 
the usual way as the mobility edge is approached from 
the localized side. We find no evidence for a diverging 
correlation length on the delocalized side. These results 
seem to agree with those of Gat and Olami but not with 
those of Galanti and Olami. We believe that our results 
are reliable, accurately portraying how vibrations look 
in our model, and perhaps generally in glasses, at ener- 
gies distant by more than 1 meV from the mobility edge. 
It is likely that for energies distant by less than 1 meV 
from the mobility edge, correlation lengths diverging with 
a universal critical exponent occur on both sides of the 
mobility edge, but this information is outside our grasp. 

The usual method for extracting information about the 
critical regime from finite systems in d = 3 is to use 
finite-size scaling combined with calculations averaged 
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over large ensembles (Kramer and MacKinnon 1993). By 
this method one is able indirectly to extract information 
from states whose localization length is comparable to 
the sample size, provided that the scaling hypothesis is 
accepted. Our method is completely different. Our sys- 
tem is large (L = 44A) and has more than 400 eigenstates 
which are well-localized within this system size. We can 
look directly at these eigenstates to see whether they de- 
cay exponentially. Averaging over groups of states in 
an energy interval of 1 meV enhances the quality of the 
measured decay constants, but even single states can be 
examined and measured. For the case of glasses with 
only a few percent of states localized, the mobility edge 
is easy to locate (to within 1 meV) without ensemble av- 
eraging. This was also found by Sheng, Zhou, and Zhang 
(1994). The direct proceedure is complementary to, and 
provides a necessary reality check on the usual scaling 
methods. It does not work well for simpler model sys- 
tems or systems with very large disorder. For example, 
Feldman et al. (1993) artificially pushed uj c deeper into 
the spectrum of amorphous Si by adding random mass 
disorder. At energies a few meV distant from the mo- 
bility edge, eigenstates were not sharply localized within 
the finite sample, and the clarity of location of lo c was 
lost. 

We study amorphous silicon in harmonic approxima- 
tion, using a realistic model which we have already 
studied extensively. Many of these studies are summa- 
rized in a paper hereafter referred to as I (Allen et al. 
1999). Amorphous Si has approximate tetrahedral co- 
ordination, and can be called an "over-constrained net- 
work glass." Even though it is not an especially typical 
glass, we suspect the vibrational properties under dis- 
cussion here have much in common with other glasses. 
Atomic coordinates were generated by Wooten using the 
algorithm of Wooten, Weiner, and Weaire (1985), with 
4096 atoms in a cubic box of side L = 44A, contin- 
ued periodically to infinity. The coordinates were re- 
laxed to a local minimum of the Stillinger- Weber po- 
tential (Stillinger and Weber, 1985), which is then used 
to find the harmonic restoring forces. Eigenvectors and 
eigenfrequenc ies were computed on the Linux cluster 
Galaxy ( ref: ht tp : / / www . galaxy, ams . sunysb . edu ) using 



the Scalapack (ref: http: / / www.netlib.org/ scalapack ) 
routine PDSYEV. The calculation used 20 Pentium II 
(400 Mhz) processors and required CPU time of 4.5 hours 
(including verification of the correctness of the result). 
The distribution of eigenfrequencies (which are all real) 
is shown in Fig. 1 of I. Finkemeier and von Niessen (1998) 
have made a systematic study of this model and several 
related ones, including more detailed size-dependence 
than we have done, as well as some variation with disor- 
der. Their results are for vibrational densities of states 
and localization properties are very similar to ours where 
they overlap. In addition, they discovered a surprising 
insensitivity of the location of the mobility edge to the 
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FIG. 1. Participation ratio for normal modes of the 
4096-atom model (higher lying data, plotted as squares) and 
for the 1000-atom model (lower lying data, plotted as cir- 
cles) plotted versus mode energy. The thick solid line seen 
at higher energy is the scaling result using the dashed line 
from Fig. 2 as discussed in the text. The inset shows the 
higher frequency data for the 4096-atom model on a loga- 
rithmic frequency scale, and a straight line whose equation is 
p oc (lu — 72meV) _3p with exponent p = 0.6. 

degree of amorphous disorder introduced in the model. 

The definition of a localized state is exponential decay 
of the eigenvector with distance from some center Rq: 



\e t {R n )\ 2 = | ai (i4)| 2 exp(-2|i4 - R \/^) 



(1) 



where a,i(R n ) isolates a factor which is not systematic 
but instead fluctuating in magnitude and sign. These 
fluctuations are random in the sense that one expects an 
ensemble of similarly prepared systems to show eigen- 
states for which the localization length ^ has "univer- 
sal" behavior near the mobility edge, and the deviations 
cii — a from exponential fall-off to be "pseudo-random" 
with zero mean (a = 0) and distribution P(a) indepen- 
dent of R n . Representative plots are shown in Fig. 5 of 
I. In order to confirm numerically that exponential decay 
is occuring, it is desirable to have \ei 
least 100 over the maximum available 
is L/2 or 22 A This means that accurately determined 
values of £ are 8A or less. Very large system sizes are 
needed to measure the large values of £ expected close to 
the mobility edge. 

Our previous studies on various size models have shown 
that most properties are reproducible from model to 
model. Figure 1 shows how reproducible is the partic- 
ipation ratio pi defined as (Bell and Dean, 1970) 



2 decrease by at 
\R n — Rq\ which 



Pi 



(2) 



The numerator is simply the normalization factor, 
equalling 1 for normalized vectors. If the mode i is a 
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pure plane wave, the denominator sums to l/N, so pi is 
the number of atoms, N. If the mode i is localized on 
a single site, the denominator is 1 and pi — 1. Roughly 
speaking, pi measures the number of atoms which par- 
ticipate in the vibration. To be more quantitative, the 
normalization integral can be written as 

1 = J drp(r)\ ai (r)\ 2 e~ 2r /Z (3) 

where the integral goes over the volume V of a 4096-atom 
cell, and the density p(f) is J2 n 5(r—R n )- Approximating 
the density as N/V and |ai(r)| 2 by its average a 2 , for 
a state decaying well within one cell volume we obtain 
a 2 w (L/£) 3 /irN. Making the same approximations to 
evaluate the denominator of Eq.(^), we obtain 

p « 8TTN(£/L) 3 a 2 ~ 2 /7? [ « {8nN/3)((/Lf, (4) 

where the assumption is made that the random numbers 
di have a Gaussian distribution. 

According to scaling theory, eigenvectors for states 
near the mobility edge (which serves as the critical point) 
should be characterized by a correlation length. On the 
localized side of the mobility edge, the decay length serves 
as the correlation length. The correlation length should 
diverge upon approaching the mobility edge. In the "crit- 
ical region" the divergence is characterized by an expo- 
nent v, 

£ oc \u> — io L \~ v => p cx \u> — oj c \~ 3i/ . (5) 

The exponent v is supposed to be "universal," that is, the 
same for all systems with time-reversal symmetry, and is 
found to have the value v = 1.57 ± 0.02 by numerical 
study of electron models in three dimensions (Ohtsuki et 
al, 1999). The data for are sufficiently dense that we 
can test this. We choose only modes with u>i >72.5meV, 
with pi <300. Making log-log plots of p versus ui — uj c for 
various choices of oj c (see the inset to Fig. 1), the data are 
roughly consistent with power laws, but precision is low 
and it is not possible even to decide where the mobility 
edge lies with much accuracy, much less to find a critical 
exponent. The data fit roughly to £ oc (u> — 69meV)~ 15 , 
(cj-70meV)- 1 ' 2 , (w-71mcV)-°- 75 , and (w-72mcV)-°- 6 . 
Therefore we try a different method which uses more of 
the information in the eigenvectors. 

For each mode in the spectrum we have located the 
atom Rq on which |ei| 2 is maximum, and computed the 
correlation function 

fi(R) = J2 |el(i?n)| 2 ?(i? - \R n - Ro\)/J2^( R ~ l#* - 

n n 

(6) 

where 5 (x) is a rectangularly shaped approximation to 
a delta function with width 0.2A The correlation func- 
tions fi(R) were then averaged over all states i with u>i 
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FIG. 2. Decay constant a = 2/£ in reciprocal Angstroms 
versus energy u) in meV for vibrations in amorphous silicon. 
Dashed lines show scaling behavior with exponent u=l. The 
curved solid line shows scaling behavior with ^=1.57 

in intervals (to, to + du>) of width du=l meV. This means 
typically 100 states were averaged. These averaged corre- 
lation functions f(uj,R) were then plotted semilogarith- 
mically over the spatial range < R < 20A. For all 
u > 72 meV, the semilog graph fell linearly (with ran- 
dom fluctuations) over a range R\ < R < 20 A, and decay 
constants a and errors Aa could be found. The distance 
Ri below which non-exponential behavior occurred was 
about 13A for the frequency range 72-73 meV, just above 
u) c , and decreased steadily to 4 A as w increased. Values 
of a and errors are shown in Fig. 2. 

The data for u) >72meV are consistent with a linear re- 
lation 2/f = 0.094(1 meV)- 1 ^ - 72 meV), that is, with 
exponent v = 1. These data are also consistent with the 
curve 2ji = 0.03(1 meV 1,57 ) _1 (a; - 71 meV) 157 , that 
is, with exponent v = 1.57. These results can be used 
in Eq.(4) to model the behavior of pi for uji > 72meV. 
The linear fit is shown as the solid line on Fig. 1. A 
much improved fit to the data in Fig. 1 is obtained by 
shifting the mobility edge down from 72meV to 70mcV. 
The disagreement is surprisingly large. Using the fit with 
v = 1.57 provides no improvement. The explanation is 
probably that pi is dominated by exponential decay (as 
assumed in Eq. (^)) only very near u c where larger sys- 
tem sizes are needed. This problem is avoided when we 
directly examine the spatial behavior of the eigenstates. 
) The correlation length for delocalized states gives the 
length scale over which fluctuations in magnitude of the 
eigenvector decay back to the average value, 

\?i(Rn)\ 2 = Mi?„) exp(-|i?„ - Mol/ti) + b t (R n )\ 2 (7) 

where ai(R n ) and bi(R n ) are both pseudorandom with 
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zero mean. For a normalized state, bf is 1/N. The same 
correlation function (Eq.(6)) was computed and averaged 
for states in the range 55-72 meV. The results agree qual- 
itatively with the form 



f(R) = a 2 exp(-2i?/0 + 1/N 



(8) 



which is expected if a, and bi are uncorrelated. The decay 
constants a = 2/£ are plotted on Fig. 2. The value 
of £ seems fixed near 7 A. This result contradicts the 
expectation of scaling theory that £ should diverge with 
the same exponent on both sides of the mobility edge. 

We believe that the numerical data of Fig. 2 represent 
reliable values of exponential decay lengths for our model 
of amorphous Si. In particular, these data show clearly 
that the mobility edge coincides with a point of apparent 
divergence of this length on the localized side, and also 
clearly show an asymmetry between the localized and de- 
localized side of lo c . Unfortunately, there is no way from 
this study to decide whether the measured lengths are 
close enough to the critical region to test the applica- 
bility of scaling theory to this system. Nevertheless, the 
surprising asymmetry seems to be a real property, similar 
to what was seen by Galanti and Olami (1999). 
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